#ifndef MLTOOL_SYMMETRIC_MATRIX_H
#define MLTOOL_SYMMETRIC_MATRIX_H

#include"layout.h"
#include"vector.h"
#include"tags.h"
#include<algorithm>
#include<functional>
#include<utility>

namespace math{
template<class T, class TRI>
class symmetric_matrix : public vector<T>
{
public:
    typedef vector<T> Array;
    typedef typename Array::difference_type difference_type;
    typedef typename Array::size_type size_type;
    typedef typename Array::value_type value_type;
    typedef typename Array::pointer pointer;
    typedef typename Array::const_pointer const_pointer;
    typedef typename Array::reference reference;
    typedef typename Array::const_reference const_reference;
    typedef typename Array::iterator iterator;
    typedef typename Array::const_iterator const_iterator;
    typedef typename Array::reverse_iterator reverse_iterator;
    typedef typename Array::const_reverse_iterator const_reverse_iterator;
    typedef matrix_tag category;
    typedef typename TRI::triangular_category triangular_category;

    template<class U, class TRI2>
    friend class symmetric_matrix;

    using vector<T>::data_;

    explicit symmetric_matrix(size_type n = 0)
        :dimension(n), vector<T>(TRI::storage_size(n, n))
    {}

    explicit symmetric_matrix(const std::pair<size_type, size_type> &dim)
        :dimension(dim.first),
          vector<T>(TRI::storage_size(dim.first , dim.first))
    {}


    symmetric_matrix(size_type n, const T &init):
        dimension(n),vector<T>(TRI::storage_size(n, n), init)
    {}



    //conversion
    template<class U>
    symmetric_matrix(const symmetric_matrix<U, TRI> & m):
        dimension(m.size1()), vector<T>(m)
    {}

    //swap
    void swap(symmetric_matrix &m)
    {
        std::swap(dimension, m.dimension);
        vector<T>::swap(m);
    }

    friend inline void swap(symmetric_matrix &m1, symmetric_matrix &m2)
    {
        m1.swap(m2);
    }

    inline T& operator()(size_type i, size_type j)
    {
        return data_[TRI::element(i, j, dimension, dimension)];
    }

    inline const T& operator()(size_type i, size_type j) const
    {
        if(TRI::one(i,j))
            return one_;
        else if(TRI::other(i,j))
            return data_[TRI::element(i,j,dimension,dimension)];
        else if(i == j)
            return zero_;
        else
            return data_[TRI::element(j,i,dimension,dimension)];
    }

    inline const size_type size1() const { return dimension; }
    inline const size_type size2() const { return dimension; }




	//operator +
    friend inline symmetric_matrix operator +(const symmetric_matrix &m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "1st called!" << std::endl;
        symmetric_matrix ans(m1.dimension);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] + m2.data_[i];
        return ans;
    }

    friend inline symmetric_matrix&& operator +(symmetric_matrix &&m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] += m2.data_[i];
        return std::move(m1);
    }

    friend inline symmetric_matrix&& operator +(const symmetric_matrix &m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] += m1.data_[i];
        return std::move(m2);
    }

    friend inline symmetric_matrix&& operator +(symmetric_matrix &&m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] += m2.data_[i];
        return std::move(m1);
    }

	//operator -
    friend inline symmetric_matrix operator -(const symmetric_matrix &m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "1st called!" << std::endl;
        symmetric_matrix ans(m1.dimension);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] - m2.data_[i];
        return ans;
    }

    friend inline symmetric_matrix&& operator -(symmetric_matrix &&m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] -= m2.data_[i];
        return std::move(m1);
    }

    friend inline symmetric_matrix&& operator -(const symmetric_matrix &m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] = m1.data_[i] - m2.data_[i];
        return std::move(m2);
    }

    friend inline symmetric_matrix&& operator -(symmetric_matrix &&m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] -= m2.data_[i];
        return std::move(m1);
    }

	//operator *
    friend inline symmetric_matrix operator *(const symmetric_matrix &m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "1st called!" << std::endl;
        symmetric_matrix ans(m1.dimension);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] * m2.data_[i];
        return ans;
    }

    friend inline symmetric_matrix&& operator *(symmetric_matrix &&m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] *= m2.data_[i];
        return std::move(m1);
    }

    friend inline symmetric_matrix&& operator *(const symmetric_matrix &m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] *= m1.data_[i];
        return std::move(m2);
    }

    friend inline symmetric_matrix&& operator *(symmetric_matrix &&m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] *= m2.data_[i];
        return std::move(m1);
    }

	//operator /
    friend inline symmetric_matrix operator /(const symmetric_matrix &m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "1st called!" << std::endl;
        symmetric_matrix ans(m1.dimension);
		for(size_type i = 0; i < m1.size(); ++i)
			ans.data_[i] = m1.data_[i] / m2.data_[i];
        return ans;
    }

    friend inline symmetric_matrix&& operator /(symmetric_matrix &&m1, const symmetric_matrix &m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "2nd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] /= m2.data_[i];
        return std::move(m1);
    }

    friend inline symmetric_matrix&& operator /(const symmetric_matrix &m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "3rd called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m2.data_[i] = m1.data_[i] / m2.data_[i];
        return std::move(m2);
    }

    friend inline symmetric_matrix&& operator /(symmetric_matrix &&m1, symmetric_matrix &&m2)
    {
#ifdef CHECK_DIMENSION_MATCH
       assert(m1.size1() == m2.size1());
       assert(m1.size2() == m2.size2());
#endif
        std::cout << "4th called!" << std::endl;
		for(size_type i = 0; i < m1.size(); ++i)
			m1.data_[i] /= m2.data_[i];
        return std::move(m1);
    }

	//operator - (负号)
	friend inline symmetric_matrix operator -(const symmetric_matrix &m)
	{
		symmetric_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = -m.data_[i];
		return ans;
	}

	friend inline symmetric_matrix&& operator -(symmetric_matrix &&m)
	{
		m *= -1;
		return std::move(m);
	}

    // 矩阵和标量的运算
    //operator +
	friend inline symmetric_matrix operator +(const symmetric_matrix &m, value_type t)
	{
		symmetric_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] + t;
		return ans;
	}

	friend inline symmetric_matrix operator +(value_type t, const symmetric_matrix &m)
	{
		return m + t;
	}

	friend inline symmetric_matrix&& operator +(symmetric_matrix &&m, value_type t)
	{
		m += t;
		return std::move(m);
	}


	friend inline symmetric_matrix&& operator +(value_type t, symmetric_matrix&&m)
	{
		m += t;
		return std::move(m);
	}

    // operator -
	friend inline symmetric_matrix operator -(const symmetric_matrix &m, value_type t)
	{
		symmetric_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] - t;
		return ans;
	}


	friend inline symmetric_matrix operator -(value_type t, const symmetric_matrix &m)
	{
		symmetric_matrix ans(m.size1(). m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = t - m.data_[i];
		return ans;
	}

	friend inline symmetric_matrix&& operator -(symmetric_matrix &&m, value_type t)

	{
		m -= t;
		return std::move(m);
	}

	friend inline symmetric_matrix&& operator -(value_type t, symmetric_matrix&&m)

	{
		for(size_type i = 0; i < m.size(); ++i)
			m.data_[i] = t - m.data_[i];
		return std::move(m);
	}
	
	// operator *
	friend inline symmetric_matrix operator *(const symmetric_matrix &m, value_type t)
	{
		symmetric_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] * t;
		return ans;
	}

	friend inline symmetric_matrix operator *(value_type t, const symmetric_matrix &m)
	{
		return m * t;
	}

	friend inline symmetric_matrix&& operator *(symmetric_matrix &&m, value_type t)
	{
		m *= t;
		return std::move(m);
	}


	friend inline symmetric_matrix&& operator *(value_type t, symmetric_matrix &&m)
	{
		m *= t;
		return std::move(m);
	}

	//operator /
	friend inline symmetric_matrix operator /(const symmetric_matrix &m, value_type t)
	{
		symmetric_matrix ans(m.size1(), m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = m.data_[i] / t;
		return ans;
	}


	friend inline symmetric_matrix operator /(value_type t, const symmetric_matrix &m)

	{
		symmetric_matrix ans(m.size1(). m.size2());
		for(size_type i = 0; i < m.size(); ++i)
			ans.data_[i] = t / m.data_[i];
		return ans;
	}


	friend inline symmetric_matrix&& operator /(symmetric_matrix &&m, value_type t)

	{
		m /= t ;
		return std::move(m);
	}

	friend inline symmetric_matrix&& operator /(value_type t, symmetric_matrix&&m)

	{
		for(size_type i = 0; i < m.size(); ++i)
			m.data_[i] = t / m.data_[i];
		return std::move(m);
	}

private:
    size_type dimension;
    static const T zero_;
    static const T one_ ;
};
template<class T, class TRI>
const T symmetric_matrix<T, TRI>::zero_ = T();

template<class T, class TRI>
const T symmetric_matrix<T, TRI>::one_ = static_cast<T>(1);

}

#endif // SYMMETRIC_H
